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Abstract.  The  structure  of  an  infinitely  strong  shock  wave  (i.e.,  a  shock  wave  with  infinitely  large  upstream 
Mach  number)  is  investigated  on  the  basis  of  the  Boltzmann  equation.  The  velocity  distribution  function 
is  expressed  as  a  sum  of  a  multiple  of  the  Dirac  delta  function,  centered  at  the  upstream  bulk  velocity, 
and  a  remainder.  Strong  evidence  that  the  remainder  has  a  singularity  in  the  molecular  velocity  space 
was  provided  by  a  previous  Monte  Carlo  simulation  for  a  hard-sphere  gas  [Cercignani  et  al.,  Phys.  Fluids 
11,  2757  (1999)].  Then,  the  singularity  was  confirmed  and  clarified  with  sufficient  accuracy  by  a  precise 
numerical  analysis  by  means  of  a  finite-difference  method.  More  specifically,  the  equation  for  the  remainder, 
which  contains  the  linear  collision  term  linearized  around  the  delta  function  and  the  nonlinear  collision 
term,  is  solved  numerically  for  a  hard-sphere  gas  after  the  nonlinear  collision  term  is  replaced  by  the  BGK 
collision  model.  The  present  paper  reports  on  the  main  result  of  this  analysis. 


INTRODUCTION 

H.  Grad  [1]  suggested  that  the  limit  of  the  shock  profile  for  the  upstream  Mach  number  going  to  infinity  exists 
(at  least  for  collision  operators  with  a  finite  collision  frequency)  and  is  given  by  a  multiple  of  the  delta  function, 
centered  at  the  upstream  bulk  velocity,  plus  a  comparatively  smooth  function,  i.e.,  the  velocity  distribution 
function  /  can  be  decomposed  as  /  =  pg6(£  —  u_)  +  /,  where  S  is  the  Dirac  delta  function,  £  is  the  molecular 
velocity,  and  u_  is  the  flow  velocity  at  upstream  infinity.  The  equation  for  the  remainder  /,  which  is  not  hard 
to  derive,  seems  to  be  more  complicated  than  the  Boltzmann  equation  itself,  but  the  presumed  smoothness 
of  its  solution  should  allow  a  simple  approximate  solution  to  be  obtained.  Grad  investigated  the  simplest 
choice  for  the  smooth  remainder,  a  Maxwellian  distribution,  the  parameters  of  which  are  determined  by  the 
conservation  equations.  Recently  the  problem  was  revisited  by  Cercignani  et  al.  [2] ,  who  provided  a  survey  on 
the  shock  wave  problem  with  particular  attention  to  an  infinitely  strong  shock  and  showed  that  a  Monte  Carlo 
simulation  for  a  hard  sphere  gas  provides  strong  evidence  for  a  singularity  of  the  remainder  /  in  velocity  space. 
The  singularity  appears  to  be  given  by  |£  —  u_|-1. 

This  singularity  was  the  object  of  a  further  investigation  by  a  deterministic,  more  accurate  method  presented 
by  Takata  et  al.  [3].  They  obtained  a  deterministic  (rather  than  Monte  Carlo)  numerical  solution  for  a  hard- 
sphere  gas,  which  confirms  the  previous  results  with  considerable  accuracy.  The  method  used  in  their  paper 
is  based  on  the  idea  that  the  singularity  is  essentially  determined  by  the  Boltzmann  collision  operator  L$ 
linearized  about  the  delta  function,  first  introduced  by  Grad  [1]  and  studied  by  Caflisch  [4] .  That  is,  Takata  et 
al.  replaced  the  nonlinear  collision  term  for  the  remainder  /  by  a  BGK  collision  model  [5,6],  keeping  the  hard- 
sphere  interaction  in  L( 5.  This  replacement,  which  simplifies  the  numerical  procedure  dramatically,  enabled 
them  to  perform  an  accurate  numerical  analysis. 

The  interest  of  the  problem  lies  in  the  fact  that,  as  conjectured  by  Cercignani  et  al.  [2] ,  the  aforementioned 
singularity  depends  on  the  molecular  model.  If  the  remainder  /  multiplied  by  the  molecular  velocity  component 
£_i_  perpendicular  to  the  direction  of  gas  flow  is  integrated  over  the  full  range  of  the  component,  and  if  the 
result  is  plotted  as  a  function  of  the  component  parallel  to  the  flow,  the  singularity  manifests  itself  as  a 
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corner  point  at  the  upstream  bulk  velocity  u_ .  The  angle  at  the  corner  point  would  change  depending  on  the 
molecular  model.  This  suggests  that  information  on  the  intermolecular  potential  might  be  obtained  from  the 
measurement  of  the  (integrated)  velocity  distribution.  For  this  reason,  we  have  started  the  research  project 
to  clarify  the  structure  of  the  singularity  for  various  molecular  models  and  obtained  the  result  for  hard-sphere 
molecules.  In  the  present  paper,  as  the  first  report  of  the  project,  we  summarize  the  result  for  hard-sphere 
molecules,  giving  a  formulation  for  more  general  molecular  models. 


BASIC  EQUATION 


We  consider  a  steady  flow  of  a  gas  through  a  standing  normal  shock  wave  and  investigate  the  structure  of 
the  shock  wave  in  the  extreme  case  of  infinitely  large  upstream  Mach  number  on  the  basis  of  the  Boltzmann 
equation.  In  this  case,  the  upstream  Maxwellian  reduces  to  p~5(£  —  u_),  where  p-  is  the  density  of  the 
gas  at  upstream  infinity.  To  analyze  this  problem,  we  first  decompose  the  velocity  distribution  function  /  as 
mentioned  in  Introduction,  i.e., 


f  =  P6(x1m-u-i)+f,  a) 

where  the  X\  axis  of  the  space  coordinates  Xi  is  taken  in  the  direction  of  the  gas  flow,  p§  is  a  function  of  X\ , 
and  i  is  the  unit  vector  in  the  X\  direction  (therefore,  u_  =  n_i).  Then  the  Boltzmann  equation  leads  to  the 
following  equations  for  /  and  ps'. 


where 


^dk  =  Ls{  ~f  )  +  Q(/~’ u~  djt  =  -K/)k=«- i  w,  (2) 

Ls (/)  =  As[2Q+(<5(£  -  u- i),  /)  -  v(5(Z,  -  u_i))/],  (3) 

Q(g,  h )  =  Q+(g,  h)  -  ^ [u(g)h  +  v(h)g],  (4) 

Q+(g,  h)  =  ±J [g(C)h(Z ')  +  g(?)h(O]B{0, \V\)d9diPd£.,  (5) 

t'ig)  =  ^j  9(Um \V\)d9d^,  (6) 

£  =  £*  -  (V  •  Q)a,  €' =  €  +  (V  •  a)a,  V  =  L  -  &  (7) 


Here,  m  is  the  mass  of  a  molecule,  a  is  a  unit  vector,  B{9,  |E|)  is  a  nonnegative  function  depending  on  the 
intermolecular  potential,  9  and  ip  are,  respectively,  the  polar  and  the  azimuthal  angle  of  a.  expressed  by  the 
polar  coordinates  with  the  polar  axis  in  the  direction  of  V,  and  the  domain  of  integration  in  Eqs.  (5)  and  (6) 
is  0  <  9  <  7r/2,  0  <  ip  <  27 r,  and  the  whole  space  of  £*.  For  hard-sphere  molecules  B  =  a2 \V\  sin0cos0  with 
a  being  the  diameter  of  a  molecule,  and  for  the  Maxwellian  molecules  B  =  (3(9)  with  (3(9 )  being  a  function  of 
9  only.  The  boundary  conditions  are  as  follows: 

/  — >  0,  ps  — >  p— ,  when  X\  — >  — oo,  (8a) 

f  -  fu+  =  <2 *£+■?/’ exp  ( -  whe"  (8b) 

where  Ft.  is  the  gas  constant  per  unit  mass,  and  p+,  T+,  and  u+i  are,  respectively,  the  density,  temperature, 
and  flow  velocity  at  downstream  infinity,  which  are  related  to  the  upstream  density  and  flow  velocity  u_. 
by  the  Rankine-Hugoniot  relations 

p+/p~  =  u-/u+  =  4,  RT+/u2_  =  3/16.  (9) 
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SOME  TRANSFORMATIONS 


In  the  present  problem,  we  can  seek  the  solution  with  cylindrical  symmetry  in  the  molecular  velocity  space, 
i.e.,  in  the  following  form: 

f  =  f(X  !,&,»?),  f?  =  (4l  +  4i)1/2-  (10) 


The  explicit  expression  of  the  linearized  collision  operator  L$  (/)  in  this  case  was  obtained  by  Caflisch  [4].  That 
is,  the  Q+(5(4  —  u_i ),/)  in  L$(f)  can  be  expressed  as 


Q+(5(£  -  u-i)J)  =  2 


/oo  /»oo 
-oo  J 0 


VA1A2 


Rb(\€  -  M-i|,  V14*  -  M-i[2  -  j|  -  u-i\- 


where 


A\  =  rjr]*  -  |4  ~  w_i|2  +  (4i  -  u_)(4i*  -  «-), 

4*  (4i*>  42*)  43*)? 

and  Rb(v,  w)  is  defined  in  terms  of  B  as 

■ffs(|Vjcos0,  |Vjsin#|)  = 


A\ 


1  1 
2  m  sin  9 


=  m*  + 14  -  u-i\2  -  (4i  -  u_)(4i*  - 
v*  =  (4i  +4i)1/2, 

[B(6,\V\)  +  B(7r/2-e,\V\)]. 


u~), 


(11) 

(12a) 

(12b) 

(13) 


Therefore,  Rb{v ,  w)  =  v  for  hard-sphere  molecules.  On  the  other  hand,  v(5(£— u_i))  in  Eq.  (3)  and  fy(/)k=«- i 
in  Eq.  (2)  are  expressed  as 


2 7T  r*/2 

-  u-i))  =  —  B{6 ,  |4  -  u-\\)dQ, 

m  Jo 

47t2  r°°  r°°  W2  ~ 

K/)k=«- i  = -  /  /  /  f(Xi,£u,r]*)B(6, 14* -u_i|)77*d0d?7*cZ4i*. 

m  J-00  Jo  Jo 


(14) 

(15) 


As  mentioned  in  Introduction,  it  is  likely  that  the  singularity  occurring  in  /  is  essentially  determined  by  the 
linearized  collision  term  In  fact,  we  have  checked  that  the  linearized  problem,  obtained  by  dropping  the 

nonlinear  term  Q(f,  /)  in  Eq.  (2),  has  a  numerical  solution  which  is  qualitatively  satisfactory  (in  particular,  it 
exhibits  the  expected  singularity)  but  fails  to  reach  a  Maxwellian  downstream  for  particles  traveling  downstream 
(for  particles  traveling  upstream  this  property  is  part  of  the  boundary  conditions).  To  simplify  the  numerical 
procedure,  therefore,  Takata  et  al.  [3]  replaced  the  nonlinear  term  by  a  BGK  collision  model  [5,6]  with  a 
suitably  adjusted  collision  frequency,  keeping  the  hard-sphere  interaction  in  L5.  In  fact,  the  main  role  of  the 
quadratic  collision  term  Q(f,  /)  is  perceived  to  be  that  of  obliging  /  to  go  to  a  Maxwellian  distribution  in  the 
entire  velocity  space.  The  replacement  mentioned  above  is  as  follows: 


where 


Q(fJ)  =  Acp(fe~f ), 


fe  = 


(2nRT)3/2 


exp 


4  —  u_i[2\ 

2  RT  )’ 


(16) 


(17a) 


/OO  poo  c\_  poo  poo 

/  /(ATi , 4i , 77)77^77(^41 ,  u  =  —  /  4i/(Ai,4i,»7)f?cM4 

-00  Jo  P  J—oqJo 


T  = 


r 


00  roo 


2ir 


m  J-o o  Jo 


14  -  u_i|2/(Xi, 4i, 77)77^77^41- 


(17b) 


Here,  Ac  is  a  constant  chosen  in  such  a  way  that  Acp+  is  the  mean  collision  frequency  of  the  molecules  in 
the  equilibrium  state  at  rest  with  density  p+  and  temperature  T+  [Ac  =  4  ttRT+  (a2 /m)  for  hard-sphere 
molecules] . 
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RESULT  OF  NUMERICAL  ANALYSIS 


The  boundary- value  problem  to  be  solved  consists  of  Eqs.  (2)  and  (3),  with  the  expressions  (11),  (14),  and 
(15)  and  with  the  replacement  (16),  and  the  boundary  conditions  (8a)  and  (8b).  It  was  analyzed  numerically 
for  hard-sphere  molecules  by  means  of  an  accurate  finite-difference  method  in  [3] .  We  summarize  the  result  in 
the  following  subsections,  where  the  origin  X\  =  0  is  set  at  the  point  at  which  the  density  is  equal  to  the  mean 
of  the  upstream  and  downstream  densities  \p  =  ( p _  +  p+)/ 2]. 

Macroscopic  quantities 

We  first  show  the  results  for  the  macroscopic  quantities.  The  profiles  of  the  density  p,  the  flow  velocity  in 
the  X\  direction  u,  and  the  temperature  T  are  shown  in  Fig.  1,  where  p*,  u*,  and  T*  are  the  dimensionless 
quantities  defined  by 


P*  - 


P~  P~ 
P+-P-  ’ 


u  —  u+ 

U*  = - , 

U-  —  u+ 


(18) 


and  l-  is  the  mean  free  path  of  the  gas  molecules  in  the  equilibrium  state  at  rest  with  density  p_,  i.e. , 
l-  =  (\[2'ko1  p_ / m)-1 .  The  figure  also  contains  the  profiles  of  the  density  p$  of  the  delta-function  part,  the 
parallel  temperature  Tj|,  and  the  perpendicular  temperature  T±  [T  =  (T||  +2Tj_)/3;  see  Eq.  (20)  below].  More 
precisely,  p$ *,  T||*,  and  Tj_*  in  Fig.  1  are  the  following  dimensionless  quantities 


(19) 


Ty  and  Tj_  being  defined  by 


T\\  =  j~pJ (6  '  u)2fdt  T±  =  zkJ +  e3)fdC  (20) 

Figure  2  shows  the  comparison  of  the  profiles  of  p,  u,  and  T  with  those  obtained  by  the  Monte  Carlo  simulation 
in  [2],  The  transition  shown  in  Fig.  1  is  monotonic  with  respect  to  X\  for  p,  u,  and  p$,  whereas  that  for  T, 
which  is  the  sum  of  monotonic  2T_l/3  and  highly  nonmonotonic  Tj|/3,  exhibits  a  very  slight  overshoot  (0.28% 
of  T+).  This  overshoot  is  smaller  than  that  observed  in  the  Monte  Carlo  simulation  (1%  of  T+).  The  overshoot 
of  almost  the  same  magnitude  was  obtained  for  several  different  lattice  systems  as  well  as  for  some  different 
choices  of  Ac  in  the  BGK  model.  However,  we  cannot  claim  the  presence  of  the  overshoot  definitely  because 


FIGURE  1.  The  profiles  of  the  macroscopic  quantities  in  the  shock  wave.  The  p*,  u*,  and  T«  are  defined  by  Eq.  (18), 
while  pgt,  Tn,,  and  Tj_*  are  defined  by  Eq.  (19). 
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FIGURE  2.  Comparison  of  the  profiles  of  the  density,  flow  velocity,  and  temperature  with  those  by  the  DSMC 
computation.  The  p*,  u»,  and  T*  are  defined  by  Eq.  (18).  The  solid  lines  indicate  the  deterministic  numerical  result  [3], 
while  the  symbols  •,  ♦,  and  ■  indicate  the  DSMC  result  [2]. 


its  magnitude  is  close  to  the  estimated  error  in  the  computation.  Because  of  the  use  of  the  BGK  model  for  the 
nonlinear  collision  term,  the  present  result  should  not  be  expected  to  give  the  precise  profiles  for  a  hard-sphere 
gas.  Nevertheless,  it  shows  reasonable  agreement  with  the  result  of  Monte  Carlo  simulation. 

Velocity  distribution  function 

We  next  show  the  behavior  of  the  non-delta-function  part  /  of  the  velocity  distribution  function.  The  contour 
lines  (u2_/ p~)r)f  =  const  at  various  points  in  the  shock  layer  are  plotted  in  the  plane  of  £i /u_  and  tj/u _  in 
Fig.  3,  where  the  corresponding  contour  lines  of  the  downstream  Maxwellian  fu+  are  also  shown  by  dashed 
lines.  The  figure  exhibits  the  concentration  of  the  contour  lines  at  the  point  corresponding  to  the  upstream 
bulk  velocity  =  u_,  r)  =  0.  The  concentration  of  the  contour  lines  means  that  the  limiting  value  of  r/f  at 

=  U-,  r]  =  0  is  a  function  of  the  angle  9 ^  =  arctan(r7/(^i  —  u_))  (0  <  9%  <  %),  i.e., 

£  lim  V)  =  F(XU  %).  (21) 

Ci  —*u- , 

The  numerical  result  of  F(X±,  0$)  /  Fmax(Xi)  at  several  points  in  the  shock  layer  is  shown  as  a  function  of  9 5 
by  solid  lines  in  Fig.  4(a),  where  -Fmax(Xi)  =  maxo<04<,r  F(Xi,9^)  and  is  given  in  Fig.  4(b).  The  fact  that 
F(X  1, 8 ^)  turns  out  to  be  finite  for  any  8 £  proves  that  /  diverges  as  |£  —  u_i|_1.  The  fact  that  the  8 £  dependence 
is  not  given  by  a  constant  times  sin#£  [as  seen  from  Fig.  4(a)]  shows  that  the  dependence  on  the  Cartesian 
components  of  the  molecular  velocity  exhibits  a  more  complex  singularity  than  a  simple  product  of  |£  —  u_i|-1 
by  a  regular  function  of  these  components.  The  concentration  of  the  contour  lines  has  been  predicted  by  the 
Monte  Carlo  simulation  in  [2].  It  was  confirmed  in  a  clearer  way  by  the  present  deterministic  computation. 

Another  aspect  of  the  singularity  can  be  viewed  if  one  considers  the  following  marginal  velocity  distribution 
function  associated  with  /: 

/OO  pOO  poo 

/  /(V1,C)^2^3  =  27r  /  rifiX^HurOdri.  (22) 

-co  J — co  J 0 

The  (u-/p^)h  versus  £i/«-  at  various  points  in  the  gas  is  shown  in  Fig.  5.  The  singularity  discussed  above 
manifests  itself  as  a  corner  point  at  £1  =  u_.  The  angle  of  the  corner  depends  on  the  position.  It  is  relatively 
small  at  X\/l _  =  —1  and  —0.6  and  increases  as  the  position  moves  downstream.  The  corner  becomes  invisible 
and  h  approaches  the  corresponding  marginal  distribution  of  /m+  as  X\/F  becomes  large  (Xj /C  =  1.4  — >  5). 
The  presence  of  the  corner  point  at  £1  =  U-  was  also  anticipated  by  the  Monte  Carlo  simulation  (see  Fig.  2 
of  [2]).  But  because  of  the  noise  or  fluctuation  intrinsic  in  the  Monte  Carlo  simulation,  it  was  not  possible 
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(b  Xy/l 


(f)  Xxl 


(c)  Xxl 


g  X1  /l 


Xi/l-  =  -0.2 


Xx/l-  =  3 


FIGURE  3.  The  contour  lines  (it_/p_)r//(Xi,£i,7/)  =  const  in  the  £i7y-plane  for  various  positions  in  the  gas.  (a) 
Xi/l-  =  -1,  (b)  Xi /|_  =  —0.6,  (c)  Xi/l-  =  -0.4,  (d)  X1/l-  =  -0.2,  (e)  Xi/U  =  0,  (f)  Xi /!_  =  0.2,  (g) 
Xi/l-  =  1,  (h)  Xi/l-  =  3.  Here,  the  dashed  lines  indicate  the  contour  lines  corresponding  to  the  downstream 
Maxwellian  distribution  /m+-  The  contour  lines  / p~)r]f  =  0.032m  (m  =  1,  2  ...)  are  plotted  (up  to  m  =  25  for 

the  Maxwellian).  The  outermost  line  corresponds  to  m  =  1  for  both  distributions. 


FIGURE  4.  The  functions  F(Xi,0$)  and  Fmax(Xi).  (a)  F/Fmax  for  X\/l-  =  —4,  —3,  —2,  —1,  —0.4,  0,  0.4,  1  and 
2,  (b)  Fmax.  Here,  Fmax(X i)  =  max  F(X i,0t).  In  (a),  the  solid  line  indicates  the  numerical  result  of  f/Fmal,  and 

O<0£<7T 

the  dot-dash  line  indicates  sin  0$. 


to  draw  a  definite  conclusion.  The  present  computation  not  only  confirms  the  existence  of  the  corner  point 
clearly,  but  also  gives  detailed  information  about  the  corner  angle. 


CONCLUDING  REMARKS 

We  have  presented  the  results  of  numerical  analysis  [3]  of  the  structure  of  an  infinitely  strong  shock  wave  in 
a  hard-sphere  gas.  We  first  split  the  velocity  distribution  function  into  a  multiple  of  the  delta  function  (delta 
part)  and  the  remainder  (non-delta  part)  and  solved  the  equation  for  the  latter  (coupled  with  the  density  of  the 
delta  part)  numerically  by  a  finite-difference  method.  In  the  process  of  analysis,  though  the  exact  form  of  the 
linear  collision  term  for  the  interaction  between  the  delta  and  the  non-delta  part  was  retained,  the  nonlinear 
collision  term  for  the  non-delta  part  was  replaced  by  the  BGK  collision  model.  The  result  confirms,  in  a  clearer 
form,  the  singular  behavior  of  the  non-delta  part  found  by  a  Monte  Carlo  simulation  in  a  previous  paper  [2]. 
The  information  summarized  in  the  present  paper  would  give  an  insight  required  for  a  mathematically  rigorous 
analysis  of  the  singularity. 

As  mentioned  in  Introduction,  the  behavior  of  the  singularity  for  different  molecular  models  would  give 
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FIGURE  5.  The  marginal  velocity  distribution  function  /i(Xi,£i)  defined  by  Eq.  (22)  as  a  function  of  fi  for  various 
positions  in  the  gas.  (a)  —2  <  X\jl _  <  0.2,  (b)  0.2  <  Xi  //_  <  5. 


information  of  physical  interest  and  importance  (see  also  the  discussion  at  the  end  of  [2]).  Since  the  use  of  the 
BGK  model  for  the  nonlinear  collision  term  for  the  non-delta  part  turned  out  to  be  successful,  the  extension 
of  the  analysis  to  molecular  models  other  than  hard  spheres,  such  as  the  Maxwellian  molecules,  is  not  difficult. 
This  will  be  the  object  of  our  future  study. 
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